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Abstract — The  value  memristor  devices  offer  to  the 
neuromorphic  computing  hardware  design  community  rests  on 
the  ability  to  provide  effective  device  models  that  can  enable 
large  scale  integrated  computing  architecture  application 
simulations.  Therefore,  it  is  imperative  to  develop  practical, 
functional  device  models  of  minimum  mathematical  complexity 
for  fast,  reliable,  and  accurate  computing  architecture 
technology  design  and  simulation.  To  this  end,  various  device 
models  have  been  proposed  in  the  literature  seeking  to 
characterize  the  physical  electronic  and  time  domain 
behavioral  properties  of  memristor  devices.  In  this  work,  we 
analyze  some  promising  and  practical  non-quasi-static  linear 
and  non-linear  memristor  device  models  for  neuromorphic 
circuit  design  and  computing  architecture  simulation. 

I.  Introduction 

HE  neuromorphic  computing  hardware  community  has 
been  re-energized  by  the  discovery  of  the  physical 
memristor  device  by  researchers  at  Hewlett-Packard  (HP) 
Laboratories,  in  Palo  Alto,  California,  in  2008  [1].  The 
memristor  device,  whose  name  comes  from  the  contraction 
of  “memory  resistor,”  has  been  characterized  as  the 
functional  equivalent  to  the  synapse  [1].  Leon  Chua 
theorized  the  existence  of  the  memristor  device  in  1971  as 
the  fourth  basic  circuit  element  [2].  Given  the  non-volatile 
nature  of  the  memristor  device,  applications  containing  such 
devices  lay  within  memory  and  computing  applications  [1]. 
As  mentioned,  the  memristor  device  operates  analogously  to 
the  biological  synapse  [1]— [3];  therefore,  it  represents  a  step 
forward  in  the  development  of  low  power  and  large  scale 
neuromorphic  computing  hardware  and  applications. 

In  order  to  apply  memristor  device  technology  to  large 
scale  computing  systems,  it  is  important  to  accurately  model 
and  simulate  its  time  domain  electronic  characteristic 
behavior.  Memristor  devices  exhibit  a  strong  hysteresis; 
therefore,  based  on  the  current  device  resistance  (or 
memristance)  state  or  initial  conditions,  we  must  be  able  to 
accurately  predict  its  future  electronic  behavior.  Several 
models  have  been  proposed  in  the  literature  to  describe  the 
non-quasi-static  electronic  time  domain  characteristic 
behavior  of  memristor  devices  [4],  [6],  [7].  In  this  work,  we 
present  a  memristor  modeling  simulation  analysis  and 
comparison  of  published  linear  and  non-linear  closed-form 
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dynamical  memristor  device  models.  We  believe  that  a  solid 
understanding  of  memristor  modeling  and  simulation 
methodologies  will  lead  to  accelerated  design  and 
development  of  memristor  powered  technologies  such  as 
neuromorphic  computing  hardware. 

II.  Memristor  Device  Models 
A.  Linear  Boundary  Drift  Model 

The  linear  memristor  device  model  reported  by  Hewlett- 
Packard  [  1  ]  [4]  states  that  the  effective  transport  mechanism 
in  Ti02  based  memristor  devices  is  through  the  drift  of 
vacancies  originating  within  an  oxygen  deficient  Ti02_x  layer 
[4].  The  Ti02  based  memristor  devices’  physical  quasi-static 
transport  mechanisms  have  been  recently  described  in  some 
detail  by  Pickett  et  al.  [5].  As  the  oxygen  vacancies  drift 
under  an  applied  external  electric  field,  the  stoichiometric 
Ti02  will  become  doped  with  the  ionized  vacancies. 
Treating  the  doped  (oxygen  vacancy  rich  regions)  and 
undoped  regions  of  the  device  as  a  pair  of  resistors  in  series, 
the  memristance  corresponding  to  a  given  boundary  position 
w  relative  to  the  device  length  or  thickness  D  can  be 
described  as  follows  [4] : 

M(w)  =  Ron{y)+Roff(l-  ■£),  (1) 

where  Ron  is  the  resistance  of  the  doped  region  and  R0ff  is  the 
resistance  of  the  undoped  region.  Again,  this  model 
describes  the  memristor  device  as  two  resistors  in  series, 
where  R0ff  and  Ron  are  the  maximum  and  minimum 
memristance  states  achievable  by  the  device,  respectively.  A 
graphical  schematic  representation  of  the  memristor  device 
model  is  shown  in  figure  1.  from  the  figure,  we  can  observe 
the  doped/undoped  boundary  region  interface,  a  dashed  line, 
whose  position,  w,  along  the  length  of  the  device,  D ,  will 
determine  the  effective  total  memristance  state  of  the  device. 

The  drift  velocity,  vD,  at  which  the  doped/undoped 
boundary  interface  moves  is  defined  as  [6] 


where  the  oxygen  vacancies  have  a  characteristic  drift 
mobility,  fiD,  under  any  applied  bias  voltage,  rj  indicates  the 
polarity  of  the  memristor,  where  rj  =  1  or  -1  for  a  device 
whose  doped  region  is  expanding  or  shrinking  respectively 
under  an  applied  voltage  bias.  Lor  example,  the  doped  region 
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of  the  memristor  device  in  Figure  1  is  on  the  left  side;  so  the 
memristor  has  an  q  =  1  polarity. 
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Fig.  1.  Graphical  model  representation  of  the  memristor  device  as 
two  resistors  in  series. 


Integrating  both  sides  of  (2)  gives  the  doped/undoped 
boundary  position  w  as  a  function  of  time  [6] 


again  letting  q(0)  =  0.  Substituting  (12)  into  (4),  we  obtained 
an  equation  for  memristance  as  function  of  charge  [6] 


M(q)  =  Ro 


277  AR  (p{t) 
Qo  Ro 


(13) 


Finally,  we  can  insert  (13)  into  (9)  to  solve  for  the  current 
flowing  through  the  memristor  device  [6] 


I(t) 


V(t) 


Ro 


27]  AR  cp(t) 
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(14) 


w(t)  =w0  +  ~~~  <7(0,  (3) 

letting  q(0)  =  0.  Substituting  (3)  into  (1),  we  can  solve  for 
the  device’s  memristance,  M(q ),  as  a  function  of  charge  [6] 


The  linear  boundary  drift  model  assumes  that  the  oxygen 
vacancies  are  free  to  traverse  the  entire  length  of  the 
memristor  unhindered  by  the  boundary  conditions  of  the 
device.  The  utility  of  this  model  lies  within  the  ease  of  usage 
and  closed  form  solution. 


M(q)  =  R0-V-^q(t),  (4) 

where,  after  grouping  terms,  the  parameters  R0,  Qo,  and  AR 
are  given  by  [6] 

R0  =  R„n  ®  +  Ro„  (l  -  f),  (5) 

0o  =  77-’  (6) 

Pq  Kon 

and 

AR  =  R0ff~  Ron  •  (2) 

From  Chua’s  seminal  memristance  equation  [2] 

d  (p  =  M  d  q,  (8) 

one  may  derive  essentially  Ohm’s  Law, 


B.  Non-linear  Boundary  Drift  Models 
The  linear  boundary  drift  model  reproduces  the 
characteristic  time  hysteresis  behavior  of  memristor  devices; 
however,  the  model  suffers  from  oversimplifications  of  basic 
electrodynamics.  Physically,  w  could  never  reach  a  zero 
width  length  because  it  would  mean  that  there  are  physically 
no  oxygen  vacancies  present  in  the  device  to  enable  the 
charge  transport  mechanisms.  On  the  other  hand,  the  entire 
length  of  the  device  could  potentially  become  doped  with  the 
oxygen  vacancies.  Modeling  the  doped/undoped  boundary 
drift  velocity  as  a  mass  on  a  spring,  the  drift  velocity,  vD , 
should  be  greatest  at  the  center  of  the  device  and  reduced  to 
essentially  zero  as  w  approaches  either  edge  (w  =  0  and  w  = 
D).  These  boundary  value  restrictions  can  be  implemented 
by  multiplying  an  additional  w  dependant  function  to  (2)  as 
shown  below  [6]  [7] 


M(9(0)=2^  =  t|.  (9) 

Using  (4),  we  can  rewrite  (9)  as  [6] 

KC0=  [K--^r<fC0]S-  (10) 

Then  integrating  (10)  over  time,  we  can  solve  for  the 
magnetic  flux 

cp{t)  =  R0q{t) ,  (11) 

which,  in  turn,  provides  an  equation  for  q(t)  via  its  quadratic 
solution  [6] 


=  (15) 

where  the  function  F  would  be  a  window  a  function 

with  non-zero  values  over  the  interval  (0,  D).  In  addition,  the 
function  should  have  its  highest  value  at  the  center  of  the 
device  (w  =  D/2)  and  be  zero  at  the  boundaries,  w  =  0  and 
w  =  D.  Joglekar  et  al.  [6]  proposed  the  window  function 

„6, 

where  p  is  a  positive  integer.  Figure  2  displays  a  graphical 
representation  of  the  window  function  described  by  (16)  for 
various  p  solutions  (p  =  1,5,  and  10).  From  the  figure,  we 
can  observe  that  the  maximum  Fp(w/D)  value  occurs  at  the 
center  of  the  device  and  that  zero  values  are  obtained  at  the 
two  boundaries.  Also,  by  varying  the  p  parameter,  we  can 
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(18) 


control  the  rate  of  change  of  the  function.  Lower  p  values 
correspond  to  a  lower  rate  of  change  and  vice  versa. 
Inserting  (16)  into  (15),  we  can  obtain  the  modified 
boundary  drift  velocity  equation 


M(yV(ti))=Ron{^)+Roff{  1- 


I(t,+ 1) 


ntj+i) 


(19) 


(17) 


Fig.  2.  Plot  of  non-linear  window  function  proposed  by  Joglekar  et 
al.  for/?  =  1,  5,  and  10. 

We  observe  that  (17)  reduces  to  the  linear  boundary  drift 
model  described  by  (2)  as  p  approaches  infinity  [6]. 
Equation  (17)  also  utilizes  the  r]  parameter,  which  is  used  to 
specify  the  physical  orientation  of  the  memristor  device.  As 
shown  in  Figure  1,  memristor  devices  are  asymmetric 
devices.  Therefore,  during  modeling  and  simulation  it  is 
important  to  consistently  specify  the  orientation  of  each 
device’s  wiring. 


Fig.  3.  Plot  of  non-linear  window  function  proposed  by  Biolek  et  al. 
for/?  =  1,  5,  and  10. 

The  nonlinear  boundary  drift  model  described  by  (17)  is 
more  physically  accurate  when  compared  to  the  linear 
model;  however,  in  order  to  solve  for  w  as  function  of  time, 
the  window  function  makes  (17)  challenging  to  integrate  for 
an  arbitrary  p.  Therefore,  a  time- step  numerical  solutions 
approach  was  employed  for  simulations.  The  following 
formulae  were  independently  derived  from  the  algebraic 
manipulation  of  (1),  (9),  and  (17)  as  shown  below 


,  (20) 

w(ti+ 1)  =  vD(ti+ 1)  [ti+ 1  -  ti]  +  w(ti)  ,  (21) 

q(ti+ 1)  =  &(ti+ 1)/  ,  (22) 


where  ti  in  (18)  corresponds  to  the  initial  time  step  and  tm  in 
(19)  to  (22)  the  next  integral  time  step. 

The  order  of  these  time-step  equations  brought  to  light 
another  challenge  in  the  implementation  of  (16),  specifically 
when  the  doped  region  covers  the  entire  device  length  ( w/D 
=  1).  It  then  follows  that  Fp(w/D  =  1)  =  0  for  all  p,  (16). 
Thus,  w  in  (21)  does  not  change  since  vD  =  0,  (20). 
Therefore,  w/D  =  1  once  again  for  the  next  time-step  during 
simulation.  Then,  this  loop  persists  till  the  end  of  the 
simulation  without  respect  to  the  change  in  the  direction  of 
the  current,  producing  invalid  results. 

A  new  window  function  was  proposed  by  Biolek  et  al  [7] 


where 
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(o,  if  /  <  0 


(23) 


(24) 


The  window  function  is  displayed  in  Figure  3  for  various  p 
integer  values  (p  =  1,5,  and  10).  This  window  function  does 
not  model  the  boundary  drift  velocity  as  a  mass  on  a  spring. 
Rather,  the  function  is  asymmetric  in  the  way  it  limits 
changes  in  vD.  For  example,  when  w  starts  at  0,  the  function 
equals  1 .  Then,  as  w  increases,  approaching  D ,  the  function 
approaches  0  as  shown  in  Figure  3  for  p  =  1.  Once  the 
current  reverses  direction,  the  function  immediately  switches 
to  1.  Then,  as  w  decreases  back  to  0,  the  function  also 
decreases  to  0  as  displayed  in  the  figure.  When  the  current 
reverses,  the  cycle  begins  once  again.  In  order  to  compute 
vD ,  (23)  can  be  substituted  into  (20)  without  altering  the 
other  four  equations.  One  advantage  of  Biolek’ s  window 
function  is  that  it  eliminates  convergence  issues  at  the  device 
boundaries. 


III.  Results  And  Discussion 

During  model  analysis  and  simulation,  all  memristor 
models  were  simulated  in  Matlab;  and  all  bias  voltage 
sources  were  of  the  form 

V(t)  =  v0  sin (co0  t  +  0),  (25) 

where  v0  is  the  voltage  amplitude  and  0  is  an  arbitrary  phase 
shift.  Typical  simulation  input  parameter  values  are  v0  =  1  - 
5  V  and  co0  =  10  -  106  rad/s.  We  can  calculate  the  flux 
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through  the  device  as  the  time  integral  of  the  voltage  across 
it  from  (25) 

<L>(t)  =  (— )  [cos<9  -  cos(&)tf  t+0)].  (26) 

\a)0/ 


A.  Linear  Boundary  Drift  Model 


Fig.  4.  Plots  of  I(t)  &  V(t )  (a),  w(t)  (b),  V-I  hysteresis  behavior  (c), 
and  M(t)  (d)  memristor  simulation  results  using  the  linear 
boundary  drift  model,  with  parameters  juD  =  10'14  nfV'V1,  D  =  10 
nm,  x0  =  0.2,  Ron  =  1700  Q,  r  =  100,  v0  =  1  V,  co  =  2n  rad/s,  and 
V(t)  =  sin(27i  t)  V. 

The  physical  memristor  device  is  characterized  by  the 
parameters  fiD ,  w0,  D ,  R0ff,  and  Ron.  Adjustments  to  the 
dopant  mobility  parameter  directly  correlates  to  changes  in 
the  boundary  drift  velocity  as  described  in  (2).  A  slower 
(faster)  velocity  corresponds  to  smaller  (larger)  changes  in  w 
per  cycle,  which  in  turn  decreases  (increases)  the  resistance 
value  spectrum  available  to  the  memristor.  Adjusting  w0 
also  directly  alters  the  effective  range  of  resistance  values 
available  to  the  memristor.  In  general,  a  higher  w0  produces 
wider  loops  in  the  I-V  plots.  However,  neither  juD  nor  w0  can 
be  set  to  completely  arbitrary  values;  otherwise,  imaginary 
numbers  arise  in  the  equations.  Overall,  the  model  operates 
over  the  widest  range  of  parameter  values  when  the  initially 
doped  region  is  less  than  half  the  device  length.  The 
maximum  viable  jiD  and  w0  values  are  related  to  the 
frequency  of  the  voltage  source,  where  a  high  frequency 
allows  for  larger  values  in  both  parameters.  Long  devices, 
high  D  values,  display  less  memristive  effects  than  short 
devices  because,  as  is  seen  in  (4),  memristance  falls  off  as  an 
inverse  square  function. 

The  R0ff  !Ron  resistance  values  can  be  arbitrarily  set  in 
accordance  with  their  definitions.  The  ratio  r  =  R0ff  /Ron 
should  be  greater  than  10,  though  ratios  of  r  =  100  -  2000 
are  more  commonly  used.  Increasing  r  generally  reduces  the 
I-V  curve  to  a  straight  line.  Additionally,  for  any  given  D , 
hysteresis  effects  are  most  prominent  when  AR  »  R0  [6]. 

For  the  linear  dopant  boundary  drift  velocity  model, 
typical  parameters  were  juD  =  (10"12  -  10"14  m^V'^s"1),  D  = 
(10  -  50  nm),  x0  =  (0.1  -  0.6),  Ron  =  (100  -  1000  Q),  and  r  - 
(100-2000). 


Figure  4  shows  typical  simulation  results.  Figure  4(a) 
superimposes  the  input  voltage  in  time  (thin  line)  against  the 
current  in  time  (thick  line).  From  the  plot,  it  is  apparent  that 
while  the  current  lags  the  voltage,  both  curves  have  the  same 
period.  This  shows  that  the  memristor  does  not  store  any 
charge  itself  but  is  a  totally  dissipative  circuit  element  [2]. 
Figure  4(c)  depicts  the  symmetric,  smooth  hysteresis  loop  of 
an  ideal  memristor.  Figures  4(b)  &  4(d)  show  the  variation 
in  width  w  and  memristance  over  time,  respectively.  From 
the  figures,  we  can  observe  that  when  w  is  greatest, 
memristance  is  minimum  and  vice-versa.  Both  parameters 
mirror  each  other. 

B.  Non-linear  Boundary  Drift  Models 

For  modeling  and  simulation  of  non-linear  memristor 
models,  the  optimal  time-step  values,  At  =  t{+ 1  -  tt  ,  were 
determined  to  be  between  10"2  -  10"4  sec.  Then,  we 
performed  model  simulations  using  Joglekar’s  window 
function,  (16),  [6];  and  the  results  are  shown  in  Figure  5. 
From  the  results,  we  can  observe  that  the  I-V  plots,  figures 
5(a)  and  (c),  exhibit  a  more  pointed  signature  compared  to 
the  linear  model  results  in  figures  4(a)  and  (c).  While  both 
I(t )  plots  have  the  same  period  as  their  respective  voltage 
inputs,  figures  5(a)  and  (c)  are  sharper  because  of  the  usage 
of  the  window  function.  We  also  noticed,  though  not  shown 
graphically,  that  for  high  p  integer  values,  the  non-linear 
model  behaves  as  its  linear  counterpart.  It  is  important  to 
notice  that  the  memristance  and  w  plots  remain  similar  for 
both  linear  and  non-linear  models  as  shown  in  figures  4  and 

5.  Under  certain  sets  of  parameters,  the  memristor  will 
fluctuate  for  a  few  cycles  before  it  settles  on  a  consistent 
pattern  of  behavior.  However,  an  appropriate  phase  shift 
choice  eliminates  these  initial  fluctuations  as  is  shown  in  the 
results  of  Figure  5,  where  a  phase  shift  of  0.16  rad  was 
employed.  The  window  function  also  gives  the  model  added 
robustness  in  terms  of  arbitrary  parameter  range  selection.  In 
addition,  in  terms  of  parameter  selection  and  adjustment, 
both  linear  and  non-linear  models  are  similarly  affected. 

In  terms  of  simulation  stability,  certain  non-linear  model 
simulations  cannot  be  performed  for  an  arbitrary  length  of 
time  when  employing  Joglekar’s  window  function.  This 
failure  is  caused  by  the  convergence  issue  described  in 
Section  II  B.  To  partially  remedy  this  problem  for  additional 
simulation  time,  we  could  increase  D  (up  to  around  50nm, 
maintaining  physical  dimensions).  However,  it  is  not  a 
comprehensive  solution. 

In  order  to  circumvent  the  convergence  issues  originating 
from  Joglekar’s  window  function,  we  can  employ  Biolek’s 
approach  described  by  (23)  [7].  Simulation  results 

employing  Biolek’s  window  function  are  displayed  in  Figure 

6.  From  the  figures,  we  can  observe  that  the  results  preserve 
the  highly  non-linear  device  characteristic  behavior.  In 
addition,  Biolek’s  model  is  unique  because  it  allows  for 
general  asymmetric  I-V  device  behavior  modeling,  which  is 
not  realizable  except  in  extreme  circumstances  with  the  two 
previous  models.  This  is  significant  because  published 
physical  memristor  experimental  data  [4]  [5]  exhibits 
asymmetric  characteristic  behavior. 
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Fig.  5.  Plots  of  1(f)  &  V(t )  (a),  w(t)  (b),  V-I  hysteresis  behavior  (c), 
and  M(t)  (d)  memristor  simulation  results  using  non-linear  dopant 
drift  model  and  Joglekar’s  window  function,  with  parameters  juD  =  6 A 
x  1014  mV's'1,  D  =  24  nm,  x0  =  0.6,  Ron  =  100  Q,  r  =  100,  p  =  7,  v0  = 
1  V,  oo  =  8tt  rad/s,  V(t)  =  sin(£7r  t  +  0.16)  V,  &(t)  =  (^)  [cos(0.16)  - 
cos(&r  t  +  0.16)]  Wb,  and  At  =  10‘4  sec. 


functions.  Once  robust,  compact  memristor  models  are  in 
place,  circuit  level  simulations  will  allow  for  applications  to 
neuromorphic  computing  architecture  development. 
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Fig.  6.  Plots  of  1(f)  &  V(t )  (a),  w(t)  (b),  V-I  hysteresis  behavior  (c), 
and  M(t )  (d)  memristor  simulation  results  using  non-linear  dopant 
drift  model  and  Biolek’s  window  function,  with  parameters  pD  =  4.4  x 
10'13m2-V1-s‘1,£>  =  41  nm,  jc0  =  0.11  ,Ron=  100  Q,  r=  10,  p  =  7,  v0  = 
1  V,  a  =  100  rad/s,  V(t)  =  sin(100  t  +  0.62)  V,  ®(t)  =  (0.01) 
[cos(0.62)  -  cos(lOOT)],  and  At  =  10'4  sec. 


IV.  Conclusion 

In  this  work,  we  analyzed  various  published,  dynamic 
linear  and  non-linear  memristor  device  models.  From  our 
study,  we  observed  that  the  non-linear  models  offer  closer 
dynamic  device  characteristic  representations  when 
compared  to  the  limited  physical  published  results  as 
opposed  to  the  linear  model.  The  non-linear  models, 
characterized  by  unique  window  functions,  provide  insight 
into  the  dynamics  of  memristor  devices. 

Future  work  will  include  performing  model-to-hardware 
correlations  to  physical  experimental  data  when  device 
fabrication  is  completed.  This  will  provide  an  opportunity 
for  refining  the  non-linear  memristor  models  and  window 
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